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In the present article we carry out a molecular dynamics study ol the core-softened system and 
show that the existence ol the water-like anomalies in this system depends on the trajectory in 
P — p — T space along which the behavior ol the system is studied. For example, diffusion and 
structural anomalies are visible along isotherms, but disappears along the isochores and isobars, 
while density anomaly exists along isochors. We analyze the applicability ol the Rosenleld entropy 
scaling relations to this system in the regions with the water-like anomalies. It is shown that the 
validity ol the ol Rosenleld scaling relation lor the diffusion coefficient also depends on the trajectory 
in the P — p — T space along which the kinetic coefficients and the excess entropy are calculated. 
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I. INTRODUCTION 

It is well known that some liquids (for example, wa- 
ter, silica, silicon, carbon, and phosphorus) show an 
anomalous behavior in the vicinity of their freezing 
lines [H-fl~8j|. The water phase diagrams have regions 
where a thermal expansion coefficient is negative (den- 
sity anomaly), self-diffusivity increases upon compres- 
sion (diffusion anomaly), and the structural order of 
the system decreases with increasing pressure (structural 
anomaly) 0, Q . The regions where these anomalies take 
lace form nested domains in the density-temperature 
(or pressure-temperature Q) planes: the density 
anomaly region is inside the diffusion anomaly domain, 
and both of these anomalous regions are inside a broader 
structurally anomalous region. It is reasonable to relate 
this kind of behavior to the orientational anisotropy of 
the potentials, however, a number of studies demonstrate 
waterlike anomalies in fluids that interact through spher- 
ically symmetric potentials 19l-47|. 

As it was mentioned above, three types of anomalies 
are the most discussed in literature with respect to the 
core-softened fluids - diffusion anomaly, density anomaly 
and structural anomaly. Here we briefly discuss these 
anomalies. 

If we consider a simple liquid (for, example, Lennard- 
Jones liquid), and trace the diffusion along an isotherm 
we find that the diffusion decreases under densification. 
This observation is intuitively clear - if density increases 
the free volume decreases and the particles have less free- 
dom to move. However, some substances have a region in 
density - temperature plane where diffusion grows under 
densification. This is called anomalous diffusion region 
which reflects the contradiction of this behavior with the 
free volume picture described above. This means that dif- 
fusion anomaly involves more complex mechanisms which 
will be discussed below. 

The second anomaly mentioned above is density 
anomaly. It means that density increases upon heating or 



that the thermal expansion coefficient becomes negative. 
Using the thermodynamic relation (dP/dT) v — ap/Kx, 
where ap is a thermal expansion coefficient and Kt is the 
isothermal compressibility and taking into account that 
Kt is always positive and finite for systems in equilib- 
rium not at a critical point, we conclude that density 
anomaly corresponds to minimum of the pressure depen- 
dence on temperature along an isochor. This is the most 
convenient indicator of density anomaly in computer sim- 
ulation. 

The last anomaly we discuss here is structural anomaly. 
Initially this anomaly was introduced via order param- 
eters characterizing the local order in liquid. However, 
later on the local order was also related to excess entropy 
of the liquid which is defined as the difference between 
the entropy and the ideal gas entropy at the same (p, T) 
point: S ex — S — Sid- In normal liquid excess entropy 
is monotonically decaying function of density along an 
isotherm while in anomalous liquids it demonstrates in- 
creasing in some region. This allows to define the bound- 
aries of structural anomaly at given temperature as min- 
imum and maximum of excess entropy. 

The problem of anomalous behavior of core-softened 
fluids was widely discussed in literature (see, for exam- 
ple, the recent review [l||). It was shown that for some 
systems the anomalies take place while for others does 
not. In this respect the question of criteria of anomalous 
behavior appearance remains the central one. However, 
another important point is still lacking in the literature - 
the behavior of anomalies along different thermodynamic 
trajectories. Here we call as "trajectory" a set of points 
belonging to some path in (P,p,T) space. For example, 
the set of points belonging to the same isotherm we call 
as "isothermal trajectory" or shortly isotherm. 

In our previous work we showed 65j that anomalies 



can exist along some trajectories while along others the 
liquid behaves as a simple one. Taking into account this 
result, it is interesting to study the behavior of the quan- 
tities demonstrating anomalies along the different physi- 
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cally significant trajectories (isotherms, isochors, isobars 
and adiabats). This investigation will allow to get deeper 
understanding of the relations between anomalous behav- 
ior and thermodynamic parameters of the system which 
spread light on the connection between thermodynamic, 
structural and dynamic properties of liquids. 



II. SYSTEM AND METHODS 

The simplest form of core-softened potential is the so 
called Repulsive Step Potential which is defined as fol- 
lowing: 



$(r) 




(1) 



where d is the diameter of the hard core, a is the width 
of the repulsive step, and e its height. In the low- 
temperature limit T = fceT '/e << 1 the system reduces 
to a hard-sphere systems with hard-sphere diameter a, 
whilst in the limit T >> 1 the system reduces to a hard- 
sphere model with a hard-sphere diameter d. For this rea- 
son, melting at high and low temperatures follows simply 
from the hard-sphere melting curve P — cT/a' 3 , where 
c ~ 12 and a' is the relevant hard-sphere diameter (tr and 
d, respectively) . A changeover from the low-T to high-T 
melting behavior should occur for T = 0(1). The precise 
form of the phase diagram depends on the ratio s = a/d. 
For large enough values of s one should expect to observe 
in the resulting melting curve a maximum that should 
disappear as s — > 1. The phase behavior in the crossover 
region may be very complex, as shown in 43j, |47|. 

In the present simulations we have used a smoothed 
version of the repulsive step potential (Eq. (JTJ), which 
has the form: 

3(r)=f~) + ~e(l-ttmh(k (r-a s ))) (2) 

where n — 14, ko = 10. We have considered er s = 1.35. 
Here and below we refer to this potential as to smooth 
Repulsive Shoulder System (RSS). 

In the remainder of this paper we use the dimensionless 
quantities: f = r/d, P = Pd 3 /e, V = V/Nd 3 = 1/p. As 
we will only use these reduced variables, we omit the 
tildes. 

In Refs. (4^, E3| 5 phase diagrams of repulsive-step mod- 
els were reported for a s = 1.15, 1.35, 1.55, 1.8. To deter- 
mine the phase diagram at non-zero temperature, we per- 
formed NVT MD simulations combined with free-energy 
calculations. For equilibration rescaling of the velocities 
was used, for sampling we used NVE ensemble monitor- 
ing the stability of the temperature. In all cases, periodic 
boundary conditions were used. The number of particles 
varied between 250, 500 and 864. No system-size de- 
pendence of the results was observed. The system was 



equilibrated for 5 x 10 6 MD time steps. Data were sub- 
sequently collected during 3 x 10 6 <5i where the time step 
St = 5x 10~ 5 . 

In order to map out the phase diagram of the system, 
we computed its Helmholtz free energy using thermody- 
namic integration: the free energy of the liquid phase 
was computed via thermodynamic integration from the 
dilute gas limit Q , and the free energy of the solid phase 
was computed by thermodynamic integration to an Ein- 
stein crystal [H, |4t|. In the MC simulations of solid 
phases, data were collected during 5 x 10 4 cycles after 
equilibration. To improve the statistics (and to check for 
internal consistency) the free energy of the solid was com- 
puted at many dozens of different state-points and fitted 
to multinomial function. The fitting function we used 
is a p . q T p V q , where T and V — 1/p are the temperature 
and specific volume and powers p and q are connected 
through p + q < N. The value N we used for the most 
of calculations is 5. For the low-density FCC phase N 
was taken equal to 4, since we had less data points. The 
transition points were determined using a double-tangent 
construction. 

The region where thermodynamic anomalies are ex- 
pected is situated close to the region where the system 
undergoes structural arrest. As a consequence, proper 
sampling of the phase space can be problematic. To 
overcome this problem we used parallel tempering [48j . 
Instead of simulating a single system, we consider n sys- 
tems, each running in the NVT ensemble at a different 
temperature. Systems at high temperatures go easily 
over potential barriers and systems at low temperatures 
sample the local free energy minima. The idea of parallel 
tempering is to attempt Monte Carlo temperature swaps 
between the systems at different temperatures. If the low 
and high temperatures are far apart, the probability to 
accept such a swap move is quite low. For this reason, 
we use a range of 'intermediate' temperatures. During 
a parallel tempering simulation, we generate equilibrium 
configurations for the system at all the temperatures in 
the simulation. In most cases, we used 8 temperatures 
and tried to swap them 40 times. This simulation took 
almost 24 hours running on 8 processors in parallel at 
the Joint Supercomputing Center of Russian Academy 
of Sciences. 

Fig. Q] shows the phase diagrams that we obtain from 
the free-energy calculations for cr s = 1.35. In fact, the 
phase diagrams for g ? = 1.15, 1.35, 1.55, 1.8 were already 



reported in Refs. [43j,|47[. We show these phase diagrams 
here too because they provide the "landscape" in which 
possible "water" anomalies should be positioned. 

Fig. da) shows the phase diagram of the system with 
a s = 1.35 in the p — T plane. There is a clear maximum 
in the melting curve at low densities. The phase diagram 
consists of two isostructural FCC domains corresponding 
to close packing of the small and large spheres separated 
by a sequence of structural phase transitions. This phase 
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FIG. 1: (Color online). Phase diagram of the system of par- 
ticles interacting through the potential (2) with a s = 1.35 in 
p - T and P - T planes. 



point: there are some values of the diameter ratio such 
that crystalline structures are strongly unfavorable and 
the glass phase could then be stable even at very low 
temperatures. The glassy behavior in the reentrant liq- 
uid disappears at higher temperatures. 



III. RESULTS AND DISCUSSION 
Diffusion anomaly 

The diffusion anomaly of the RSS was discussed in sev- 
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43] and 



eral our previous articles. In the Refs. 
we showed that the diffusion anomaly takes place at the 
shoulder width a\ = 1.35 while in the work [H3| the 
breakdown of the Rosenfeld scaling for this system was 
demonstrated. Further discussion of the Rosenfeld scal- 
ing was reported in the work [r35l | . The main idea of this 
paper is that the anomalous diffusion behavior present 
along low-temperature isotherms while along isochors dif- 
fusion coefficient is monotonous. As a result, Rosen- 
feld scaling is valid along isochors and high-temperature 
isotherms, however, it breaks down for the isotherms with 
anomalies. It means that appearance or not of anomalies 
at some trajectory can cause physically different behavior 
of the system. As a result, following different trajectories, 
we can or can not observe some effects. It is of particu- 
lar importance since experimental works deal mostly with 
isobars and isotherms while theoretical studies - with iso- 
chors and isotherms. Taking this into account, one can 
expect that some of the effects observed along isobars are 
not visible along isotherms and isochors which makes us 
confused while comparing experimental results with the- 
oretical predictions. Here we extend the study of anoma- 
lies to four different physically meaningful trajectories: 
isotherms, isochors, isobars and adiabats. 



diagram was discussed in detail in our previous publica- 
tion 43j,|47|. It is important to note that there is a region 



of the phase diagram where we have not found any sta- 
ble crystal phase. The results of Ref. [43| suggest that a 
glass transition occurs in this region with T g — 0.079 at 
p = 0.53. The apparent glass-transition temperature is 
above the melting point of the low-density FCC and FCT 
phases. If, indeed, no other crystalline phases are stable 
in this region, the "glassy" phase that we observe would 
be thermodynamically stable. This is rather unusual for 
one-component liquids. In simulations, glassy behavior 
is usually observed in metastable mixtures, where crys- 
tal nucleation is kinetically suppressed. One could argue 
that, in the glassy region, the present system behaves 
like a "quasi-binary" mixture of spheres with diameters 
d and a s and that the freezing-point depression is analo- 
gous to that expected in a binary system with a eutectic 



Fig. [21a) shows the diffusion coefficient of RSS system 
along a set of isotherms for u\ = 1.35. One can see that 
the diffusion coefficient demonstrates anomalous behav- 
ior for the temperatures below T = 0.25. However, if 
we look at the Fig. HJb) where the same diffusion coeffi- 
cient data are arranged along isochors we do not observe 
anomalies - the diffusion is a monotonous function of tem- 
perature along isochors. However, some of the isochors 
cross. One can see from the inset that the cross of the 
isochors corresponds to the densities between p = 0.4 and 
p = 0.6. Comparing it to the isotherms we see that this is 
the region of anomalous diffusion. It means that even if 
we do not see nonmonotonous behavior of diffusion along 
isochors we can identify the presence of anomaly from 
crossing of the isochors. However, this method seems to 
be technically more difficult since we need to measure 
many points belonging to different isochors rather then 
one isotherm. 
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FIG. 2: (Color online). Diffusion coefficient of the RSS system 
along (a) isotherms and (b) isochors. The insets show (a) the 
low temperature isotherms and (b) low temperature region of 
some isochors. 



FIG. 3: (Color online). Diffusion coefficient of the RSS system 
along isobars as a function of (a) density and (b) temperature. 
The insets show anomalous regions in the corresponding co- 
ordinates. 



Isothermal and isochoric behavior of diffusion in 
anomalous region have already been discussed in our pre- 
vious publication [65| and we give these plots here for the 
sake of completeness. Fig. [3][a) shows the diffusion co- 
efficient along a set of isobars as a function of density. 
The diffusion coefficient is again monotonous. The slope 
of the curves is always negative. However, as it can be 
seen from the inset, the slope approaches infinity at low 
diffusions at pressures P = 2.0 and P — 2.5. This cor- 
responds to densities 0.45 — 0.55 inside the anomalous 
region. One can imagine that if we lower the tempera- 
tures along these isobars we can observe change of the 
slope to positive one, however, we do not have data for 
these temperatures. 

Fig. Efb) shows the diffusion coefficient along isobars 
as a function of temperature. The situation is analogous 



to the case of isochors: the curves are monotonous, how- 
ever, they intersect at low temperatures, corresponding 
to anomalous region (see the inset in Fig.[3jb)). It means 
that if we have the diffusion coefficient along isobars we 
can identify the presence of anomalies by monitoring the 
intersections of the curves. However, by the reasons dis- 
cussed above this method is not practically convenient. 

The last physically meaningful trajectory considered 
in the present work is the adiabat. This trajectory 
is defined as constant entropy curve. The entropy is 
calculated as following. We compute excess free en- 
ergy by integrating the equation of states: ^k^T = 

7& - I^Io P P(P V kBT d P'- The excess ei ^°Py 
can be computed via S ex — . The total en- 

tropy is S = S ex + Sid, where the ideal gas entropy is 
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= lln(T) - Hp) +H ^ mkB S / ^ >2 )- The last 
term in this expression is constant and is not accounted 
in our calculations. 

The behavior of entropy itself will be discussed below. 
Here we give the diffusion coefficients along the adiabats 
(Fig.Sl(a)-(c)). 

One can see from the Fig. [4] (a)-(c) that the anomaly 
takes place along adiabats in all possible coordinates 
(look, for example, S = —4.0 adiabat). It means that in 
case of adiabatic trajectory one can identify the anoma- 
lous region monitoring any of three thermodynamic vari- 
ables (P, p, T). However, this trajectory is rather difficult 
to realize in simulation or experiment. 



Density Anomaly 

As we mentioned above density anomaly corresponds 
to appearance of a minima on isochors of the system. 
The isochors are shown in the Fig. [SJa). It is evident 
from the figure that some of the isochors do demonstrate 
minima. The location of the minimum in the p~T plane 
is shown in Fig. [5][b). One can see from this figure that 
the region of the density anomaly is located inside the 
region of the diffusion anomaly which is consistent with 
the general picture corresponding to water. 

If we turn to the shape of the isobars themselves, that 
is the dependence of temperature on density at fixed pres- 
sure we do not find any traces of anomalies there (Fig. 
[6|). Like in the case of diffusion, the curves have nega- 
tive slope which approaches zero at low temperatures and 
densities corresponding to anomalous regime. However, 
the curve remains monotonous which means that there 
is no sign of density anomaly along isobars. 

Figs. [7£a) and (b) show adiabats of the RSS in p — T 
and P—T coordinates. Interestingly, the curves are again 
monotonic, but the slope of the curves at low and high 
densities (pressures) is very different. At the same time 
the middle density (pressure) curves demonstrate a con- 
tinues change from the high slope (low-density regime) 
to low slope (high-density) regime. It allows to conclude 
that even if adiabats show monotonous behavior one can 
easily identify the location of anomalous region by the 
slope change of the curves. 



Structural Anomaly 

Structural anomaly region can be bounded by using 
the local order parameters or by excess entropy minimum 
and maximum. Here we apply the definition via excess 
entropy. 

The behavior of excess entropy is qualitatively anal- 
ogous to the behavior of diffusion coefficient. Because 
of this we briefly describe it here noting that most of 
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FIG. 4: (Color online). Diffusion coefficient of the RSS system 
along adiabats as a function of (a) temperature, (b) density, 
and (c) pressure. The insets show anomalous regions in the 
corresponding coordinates. 
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the conclusions about diffusion coefficient along different 
trajectories can be applied to excess entropy as well. 

Figs. EJa) and (b) show the excess entropy along 
isotherms and isochors. Like for the diffusion coefficient, 
excess entropy demonstrates anomalous grows in some 
density range at low temperatures. At the same time 
excess entropy is monotonous along isochors. However, 
the curves for different isochors cross which indicates the 
presence of anomaly. 

The excess entropy along isobars is monotonically de- 
creasing function of density and monotonically increas- 
ing function of temperature. However, the curves cross 
at low temperatures indicating the presence of anomalies 
as it was discussed for the case of diffusion. 

It is important to note that the range of densities of 
structural anomalies is wider then the one of diffusion 



FIG. 6: (Color online). A set of isobars of RSS. 
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FIG. 7: (Color online). Adiabats of the RSS in (a) density 
temperature and (b) pressure - temperature coordinates. 
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FIG. 8: (Color online). Excess entropy of RSS fluid along (a) 
isotherms and (b) isochors. 



FIG. 9: (Color online). Excess entropy of RSS fluid along 
isobars as a function of (a) density and (b) temperature. The 
inset in (b) shows the cross of the curves at low temperatures. 



anomaly which is consistent with the literature data for 
core-softened systems. 

To summarize this section, we emphasize that the 
anomalous behavior can be seen only along some par- 
ticular trajectories in (P, p, T) space. For example, diffu- 
sion and structural anomalies are visible along isotherms 
while density anomaly - along isochors. Other meaningful 
trajectories (isobars and adiabats) can identify the pres- 
ence of anomalies via crossing of the curves, but they do 
not allow easily identify the boundaries of the anomalous 
region. 



IV. ROSENFELD SCALING 



In 1977 Rosenfeld proposed a connection between ther 



modynamic and dynamical properties of liquids 
The main Rosenfeld's statement claims that the trans- 



port coefficients are exponential functions of the ex- 
cess entropy. In order to write the exponential rela- 
tions Rosenfeld introduced reduction of the transport 
coefficients by some macroscopic parameters of the sys- 
tem 
D* 



For the case of diffusion coefficient one writes: 
D -rr—^TT— tt7t, where m is the mass of the parti- 



(fc B T/m) 1 /2 

cles. The Rosenfeld scaling rule can be written as: 



D* = A- 



(3) 



where A and B are constants. 

In his original works Rosenfeld considered hard 
spheres, soft spheres, Lennard- Jones system and one- 
component plasma [H(J |5l|. After that the excess en- 
tropy scaling was applied to many different systems in- 
cluding core-softened liquids 
" HBl 56|, binary mixtures 57 



IMlliUHIl, lkmid metals 



58|. ionic liquids 5t| 6 



network- forming liquids 52,5^, water 61 1 , chain fluids 



[62j and bounded potentials 
In our recent publication 




we showed that for 



the case of the core-softened fluids the applicability of 
Rosenfcld relation depends on the trajectory. In partic- 
ular, Rosenfeld relation is applicable along isochors, but 
it is not applicable along isotherms. 

The breakdown of the Rosenfcld relation along 
isotherms can be seen from the following speculation. 
The regions of different anomalies do not coincide with 
each other. In particular, in the case of core-softened 
fluids the diffusion anomaly region is located inside the 
structural anomaly one. It means that there are some re- 
gions where the diffusion is still normal while the excess 
entropy is already anomalous. But this kind of behavior 
can not be consistent with the Rosenfeld scaling law. 

However, from this speculation it follows that the 
Rosenfeld scaling should hold true along the trajecto- 
ries which do not contain anomalies, i.e. isochors and 
isobars. In our recent publication [r35| we considered the 
Rosenfeld relation along isotherms and isochors. Here we 
bring these trajectories for the sake of completeness and 
add the verification of the Rosenfeld relation along iso- 
bars (Fig. [TCTa) - (c)). One can see that the Rosenfeld 
relation does break down along isotherms which is con- 
sistent with the speculation above. At the same time it 
holds true along both isochors and isobars which is con- 
sistent with the monotonous behavior of both diffusion 
coefficient and excess entropy along these trajectories. 

This observation makes evident that Rosenfeld relation 
is valid only along the trajectories without anomalous 
behavior. 
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V. CONCLUSIONS 

In conclusion, in the present article we carry out a 
molecular dynamics study of the core-softened system 
(RSS) and show that the existence of the water-like 
anomalies in this system depends on the trajectory in 
P — p — T space along which the behavior of the system 
is studied. For example, diffusion and structural anoma- 
lies are visible along isotherms, but disappears along the 
isochores and isobars, while density anomaly exists along 
isochors. We analyze the applicability of the Rosenfeld 
entropy scaling relations to this system in the regions 
with the water-like anomalies. It is shown that the va- 
lidity of the of Rosenfeld scaling relation for the diffusion 
coefficient also depends on the trajectory in the P — p — T 
space along which the kinetic coefficients and the excess 
entropy are calculated. In particular, it is valid along 
isochors and isobars, but it breaks down along isotherms. 
The breakdown of the Rosenfeld relation along isotherms 
can be explained in the following way. The regions of 
different anomalies do not coincide with each other. In 
particular, in the case of core-softened fluids the diffusion 
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FIG. 10: (Color online). Rosenfeld relation for RSS along 
isochors, (b) isotherms, and (c) isobars. 
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anomaly region is located inside the structural anomaly 
one. It means that there are some regions where the dif- 
fusion is still normal while the excess entropy is already 
anomalous. But this kind of behavior can not be consis- 
tent with the Rosenfeld scaling law. 
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